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Abstract 

We present molecular dynamics simulation measurements of the viscosities of 
the Gay-Berne phenomenological model of liquid crystals in the nematic and 
isotropic phases. The temperature-dependence of the rotational and shear 
viscosities, including the nonmonotonic behavior of one shear viscosity, are in 
good agreement with experimental data. The bulk viscosities are significantly 
larger than the shear viscosities, again in agreement with experiment. 
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Ever increasing computer power has made simulations of simple, yet realistic molecular 
models of liquid crystals a feasible enterprise. With the aid of computers it is possible to 
study the effects of molecular shapes, sizes and interactions on macroscopic behavior. Three 
types of liquid crystal models have been used in simulation work. The first is the Lebwohl- 
Lasher model |IJ, a lattice model for rotators. This model can be used to study the isotropic- 
nematic transition as a rotational order- disorder transition in an effective crystalline solid. 
Another class of models which has received much attention and includes the translational 
degrees of freedom of a liquid crystal uses hard particles of various shapes which interact 
solely by excluded volume effects These hard body models exhibit very rich phase 

diagrams including smectic, columnar and cubatic phases. Finally, in recent years there 
has been considerable numerical study of the Gay-Berne (GB) system Q which is a fluid of 
point-like particles, each carrying a unit vector u which mimics the long molecular axis. The 
particles interact via an anisotropic Lennard- Jones potential which depends on the relative 
orientation and location of a pair of molecules. This system displays rich behavior like the 
hard body models, but also includes attractive forces. These latter forces play an especially 
important role in the formation of smectic phases, which are more readily formed in the GB 
system than in the hard body models.. The GB model has also been extended to include 
chirality f| and to model discotics f|. 

Previous studies of the GB system have focused on the phase diagram |7|-|§ and single- 
particle translational and rotational dynamics [|1(J. In this Letter we report on numerical 



measurements of the viscosities of the GB system. We find a number of results which indicate 
that the GB model exhibits the principal dynamical features of a real nematic liquid crystal. 
The temperature dependence of the shear viscosities and the rotational viscosity below 
the isotropic-nematic transition is qualitatively similar to what is observed in experiment 



IT] |nj , including the nonmonotonic behavior of one of the shear viscosities. The two bulk 
viscosities are an order of magnitude larger than the shear viscosities in agreement with 
ultrasonic measurements JTH] . We also find that long-time correlations of the director exhibit 



the expected Brownian motion due primarily to the finite-size of our system. This motion 



yields exceptionally good statistics for the rotational viscosity. 

The GB potential is modeled to give the best fit to the pair potential for a molecula 
consisting of a linear array of four equidistant Lennard- Jones centers with separation of 
2a Q between the first and fourth sites (subsequent work || examined different fits) The GB 
potential is given by, 
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where Ui,u 2 are unit vectors giving the orientations of the two molecules separated by the 
position vector r. The parameters e(ui,u 2 , r) and cr(ui,u 2 ,r) are orientation dependent 
and give the well depth and intermolecular separation where U = respectively. The well 
depth is written as, 

e(ui, u 2 , r) = e E v (vi 1 , u 2 )e' /1 (ui, u 2 , r) (2) 

where 

e(u 1 ,u 2 ) = (l- X 2 (u 1 -u 2 ) 2 )- 1 / 2 (3) 

and 
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The shape anisotropy parameter \ is given by 



2 I l + x'K-u,) 
I (?-m-g-n2) 2 l (4) 



X = {(^ s ) 2 -l}/{(-f + l} (5) 



where a e and a s are the separation of end-to-end and side-by-side molecules respectively. 
The parameter x' is given by 

X' = {1 - {e e /e s f^}/{\ + (e e /e s )^} (6) 



The ratio of the well depths for end-to-end and side-by-side configurations is e e /e s . 

We simulated a system of iV = 256 particles using the molecular dynamics technique. 
The ratio o e jo s was set equal to 3, the ratio e e /e s = 5 and exponents v — 1 and fi = 2 
. These values yield the best fit to the linear array of four Lennard- Jones centers [|J, 
though other values have been used P|J. We used cubic periodic boundary conditions, 
and cut off and smoothed the potential at 3.8ero- The equations of motion were solved 
using the leap-frog algorithm with integration time-step At* = 0.001 in reduced units ( 
At* = At(ma%/e y 1/2 , where m is the mass of the molecule ). To mimic the behavior of 
p-azoxyanisole (PAA) we chose Co = 2.7 A, e /kB = 400° K, and m = 4.29 x 10~ 23 g, yielding 
At = 2.4 x 10 -12 sec. Translational and rotational temperatures were controlled using 



the Nose- Hoover thermostat ||16|| . The initial configuration was generated by locating nearly 
parallel GB molecules on the sites of an fee lattice at low scaled density, p*{= Na^/V) = 0.10. 
We carried out a long run (20,000 iterations) which reduced the nematic order parameter to 
a value of 0.12, and disordered the system translationally. The system was then gradually 
compressed at constant scaled temperature T*(= k B T/e ) = 3.0 to a density p* = 0.32. 
Finally, at constant density the temperature was lowered in small steps to the isotropic- 
nematic transition temperature and below. At this stage longer simulation runs (180,000 - 
300,000 iterations) were performed to obtain the values of the viscosities. 

The dynamical description of a compressible nematic requires six viscosities: three shear 
viscosities, ^1,^2, and z/ 3 ; two bulk viscosities, z/ 4 — v 2 , and and a director rotational 
viscosity, 71. To calculate these viscosities from correlation functions of the stress tensor 
and the director, we transformed to a coordinate system where the 3-axis is parallel to the 
average orientation of the director, and the 1 and 2 axes are perpendicular to the director. 
The elements of the stress tensor are defined by: 

= + E E '•:;'/?)• «, = 1, 2, 3, (?) 

i i j>i 

where V is the volume of the system, p l is the linear momentum of molecule i, and r 4 - 7 and 
iij are respectively the relative position vector and force between molecules i and j. The 



five viscosities v\, 1/2, ^3, ^4, and associated with shear and compression are then given in 



terms of Kubo-like formulas JDJ : 
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where 



v{t) = \{an{t)+(J22{t)) (13) 

The time-correlation functions appearing in Eqs. (8)-(12) were evaluated by averaging 
over successive time origins [|l8j. The components of the stress tensor were evaluated first 
in the cubic coordinate system (with axes parallel to the sides of the cell) where we defined 
periodic boundary conditions and then transformed to the director coordinate frame. 

In order to compare our simulation results with experimental measurements of shear 
viscosities, we applied a magnetic field H along the long diagonal of simulation cube (the 
direction along which the director spontaneously orders), with XaH 2 = 1.0 in our reduced 
units, (xa is the anisotropy in the magnetic susceptibility). Experiments typically measure 
Miesowicz viscosities [ IT |: 

Vi =^3 + t7i(1- A) 2 + A71 (14) 
V2 = ^3 + ^7i(l -A) 2 (15) 



V3 = "2 (16) 



The parameter A is a reactive coefficient ]T7|] which determines the response of the director to 



shear flow. Because we have not yet performed a direct shear flow experiment we do not have 
a value of A for the GB system; in PAA |T7| A = 1.15±0.10. We have assumed that A is unity 



in the GB system and computed Miesowicz viscosities from our measurements; a 10% change 
in this value of A will not alter the qualitative features of our results. Results for the three 
Miesowicz shear viscosities are shown in Figs. |l| and 0; results for 71 alone are given in Fig. |j. 
The large value of H needed to stabilize the director motion smears the nematic-isotropic 
transition. At sufficiently high temperatures where the nematic order parameter is low, (see 
the inset of Fig. p, the system is essentially isotropic with large director fluctuations and we 
have plotted the average of the two shear viscosities 772 and 7/3 as a single isotropic viscosity. 
Below the temperature T ~ 2.0 where the nematic order parameter becomes substantial we 
note that 772 decreases at first as the temperature is lowered and then rises. This effect has 
been observed experimentally in PAA [fTT] , [m , N-(4'-methoxybenzilidene)-4-(n butyl) aniline 
(MBBA) [0, and in the cyanobiphenyl homologues |TJ]]. The viscosity r] 2 is associated 



with shear flow parallel to the director, so its value drops when the nematic order becomes 
appreciable. Its subsequent rise is not fully understood. Likewise the fact that 773 > 772 below 
the transition has also been observed [|TTl|l2ll and deep in the nematic phase the ratio 773/772 
has been found to be approximately 1.5, in both PAA and MBBA, in good agreement with 
our results. The temperature dependence of V\ or 771 has not been measured, but the value of 



Ui deep in the nematic phase of PAA is known |Tj| and is comparable to the value of 773. We 
find V\ = 1.52 and 1.044 centipoise at T* = 1.2 and 1.4 respectively. We have also measured 
the bulk viscosities at these latter temperatures and find z/4 = 19.6 and 9. 54;^ = 25.9 and 
10.74 centipoise respectively. Experimental measurements of the bulk viscosities deep in the 
nematic phase of PAA and para-azoxyphenetole (PAP) [15j also show values which are an 
order of magnitude larger than the largest shear viscosity. The temperature dependence of 
the bulk viscosities has not been explored experimentally. Numerically determining the bulk 
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viscosities is more time consuming than determining the shear viscosities, so we have not 
explored their full temperature dependence. 

The magnitudes of our viscosities are of the correct order of magnitude for what has 
been measured in PAA. The measured values of the viscosities for MBBA are consistently 
an order of magnitude higher than the corresponding quantities in PAA. Simply changing 
the values of eo, &o, and the mass in the GB potential to values appropriate to MBBA is not 
sufficient to account for this difference. Nevertheless, the temperature dependence of the 
shear viscosities is qualitatively similar in PAA and MBBA, so the present parameterization 
of the GB potential captures the essential physics. 

Correlations in the director rotational motion and the value of the rotational viscosity 
7i are best studied by measuring the correlation function: 

C(t) =< (5n(t) - 5n(0)) 2 >, (17) 

where 

8n(t) = ^( ni (t)+n 2 (t)), (18) 

and, 

= AfE a = 1,2. (19) 

i 

If the director motion is diffusive (Brownian motion), then C(t) at long times will follow: 

C(t) ~ (20) 

and the value of 71 can be extracted with excellent statistics. This behavior was observed 
in our system (see Fig. ^ at sufficiently long times, and is sensible given the magnetic field 
and especially the small system size. At shorter times the director motion obeys the form 
C(t) ~ t 2 approximately, which is to be expected for ballistic behavior. For much larger sys- 
tem sizes and small magnetic fields we would expect this latter behavior to persist to larger 



times due to the presence of massless Goldstone modes jl9|| , even though the microscopic 



director motion is no longer ballistic. This Goldstone-like behavior has been reported by 



Zhang et al. ||20|| , in Monte Carlo studies of the Lebwohl-Lasher lattice model where very 
large system sizes are readily studied. The temperature dependence of 71 is given in Fig. |] 
in the form of an Arrhenius plot, i.e., 71 oc exp(E/ksT), where E is an activation energy. 
This Arrhenius behavior has been observed experimentally in MBBA and in a mixture of 
the two isomers of p-methoxy-p'-butylazoxybenzene (N4) ||21|| . The slope of the line give the 
activation energy E ~ O.leV. This value has not been measured in PAA, but in MBBA it 
is approximately 5 times larger. As discussed above, our numerical results are consistently 
an order of magnitude lower in value than those measured in MBBA. The value of 71 has 
been measured in PAA [2^] deep in the nematic phase (T = 122°) and the result,7i ~ 6.8 



centipoise is of the same order of magnitude as our results. 
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FIGURES 

FIG. 1. The isotropic shear and two of the Miesowicz viscosities as functions of the scaled 
temperature T*, in the presence of a magnetic field. The inset shows the nematic order parameter 
S as a function of T* for the same value of magnetic field. The parameters of the GB potential 
were chosen to mimic PAA; see the text. 

FIG. 2. The third Miesowicz viscosity, r/i, which is two to three times larger than the other 
two viscosities; compare Fig. [j] 

FIG. 3. The director correlation function C(t), Eq. (|i~7|), as a function of dimensionless time, 
at T* = 1.4. The short-time behavior is fit to a straight line with slope 1.92, while the slope at 
long times is 1.00. The error bars are imperceptible on this scale. 

FIG. 4. The rotational viscosity 71 as a function of scaled temperature indicating the activated 
form of the dependence. 
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